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A general approach to calculate the diabatic surfaces for electron-transfer reactions is presented, 
based on first-principles molecular dynamics of the active centers and their surrounding medium. 
The excitation energy corresponding to the transfer of an electron at any given ionic configuration 
(the Marcus energy gap) is accurately assessed within ground-state density-functional theory via a 
novel penalty functional for oxidation-reduction reactions that appropriately acts on the electronic 
degrees of freedom alone. The self-interaction error intrinsic to common exchange-correlation func- 
tionals is also corrected by the same penalty functional. The diabatic free-energy surfaces are then 
constructed from umbrella sampling on large ensembles of configurations. As a paradigmatic case 
study, the self-exchange reaction between ferrous and ferric ions in water is studied in detail. 
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A wide variety of processes and reactions in electro- 
chemistry, molecular electronics, and biochemistry have 
a common denominator: they involve a diabatic elec- 
tron transfer process from a donor to an acceptor 
These reactions cover processes and applications as di- 
verse as solar-energy conversion in the early steps of 
photosynthesis, oxidation-reduction reactions between a 
metallic electrode and solvated ions, and the TV char- 
acteristics of molecular-electronics devices 0. The key 
quantities of interest are the reaction rates (or, equiva- 
lently, the conductance) and the reaction pathways. Re- 
action rates, in the general scenario of Marcus theory 
0- Di 01 ■ have a thermodynamic contribution (the classi- 
cal Franck-Condon factor, broadly related to the free en- 
ergy cost of a nuclear fluctuation that makes the donor 
and the acceptor levels degenerate in energy), and an 
electronic-structure, tunneling contribution (the Landau- 
Zener term, related to the overlap of the initial and final 
states) . 

We argue in the following that state-of-the-art first- 
principles molecular dynamics calculations, together with 
several algorithmic and conceptual advances, are able to 
describe with quantitative accuracy these diabatic pro- 
cesses, while including the realistic description of the 
complex environment encountered, e.g. in nanoscale de- 
vices or at the interface between molecules and metals. 

Fig. H shows schematically an electron-transfer process 
and the free-energy diabatic surfaces according to the 
picture that was pioneered by Marcus 0, S IE ■ I n a 
polar solvent, the electron transfer process is mediated 
by thermal fluctuations of the solvent molecules. In the 
reactant state, the transferring electron is trapped at the 
donor site by solvent polarization; transfer might occur 
when the electron donor and acceptor sites become de- 
generate due to the thermal fluctuations of the solvent 
molecules. To characterize the role of the solvent on the 
electron-transfer reaction, a reaction coordinate e for a 
given ionic configuration is introduced, as the energy dif- 
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FIG. 1: Diabatic free-energy surfaces for ferrous-ferric elec- 
tron transfer reactions. AG is the free energy barrier, and A 
is the reorganization energy. The reaction coordinate e is the 
Marcus energy gap. 



ference between the product and reactant state at that 
configuration This definition of reaction coordinate 
captures the collective contributions from the solvent. 

There have been a number of pioneering classical 
molecular dynamic studies 0, |jj 0] studying the reac- 
tions between aqueous metal ions. However, quantita- 
tive agreement has not been achieved with classical force 
fields. The reorganization energy A (i.e. the free energy 
cost to reorganize the solvent molecules from the config- 
urations at equilibrium with the product to the configu- 
rations at equilibrium with the reactant without electron 
transfer) for the aqueous Fe 2+ -Fe 3+ self-exchange reac- 
tion was found to be 3.6 eV for ions 5.5 A apart (lcj . 
while experimentally [rH ] (at the slightly shorter separa- 
tion of 5.32 A) it is found to be 2.1 eV. Although there 
have been studies of electron-transfer reactions includ- 
ing electronic polarization in classical force-field poten- 
tials 0, El , full first-principles studies are required to 
describe realistically and quantitatively these reactions. 
Recently, an elegant grand-canonical density functional 
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approach has been introduced to address this class of 
problems This approach is, however, targeted 

at half-reactions for a donor or an acceptor in contact 
with an electron reservoir. 

In this paper, we present a novel technique to study 
electron-transfer reactions from first-principles molecular 
dynamics, with ferrous- ferric self-exchange as a paradig- 
matic example. We use Car-Parrinello molecular dynam- 
ics [H E3 andspin-polarized DFT in the PBE-GGA 
approximation |lj| • Fig. [21 shows schematically the sam- 
pling procedure used, following the lines of Ref. ^(| f° r 
classical simulations. An ionic trajectory is first gener- 
ated with the ions in the (2+r) and (3-r) states of charge, 
respectively, r is an "umbrella-sampling" parameter used 
to explore different regions of the phase space. We then 
perform two separate runs with the electronic state con- 
strained in the reactant or in the product configuration, 
and with the ions following the afore-generated ionic tra- 
jectories. The reaction coordinate e at every time step is 
thus given by the difference between the energies of the 
product and the reactant state. The probability distri- 
bution P(e) is then calculated as: 



Sampling run 



P(e) 



E^(r),eexp{-P[E r (r) - E s (r)]} 

T 

J2exp{~p[E r (T)~E s (T)]} ' 



(1) 



where e'(r) = E p (t) — E r (r) is the reaction coordinate 
at time r; E r (r), E p (t) and E s (t) are the energies of 
the system in the reactant, product and sampling oxi- 
dation states, respectively, at time r, and ^ e '( T ),e is the 
Kronecker delta. The exponential term in the expres- 
sion above restores the correct thermodynamical sam- 
pling according to the energy surface E r . The free energy 
G(e) is derived from the probability distribution P(e) as 
G(e) = -fc B Tln(P(e)). 

It is of central importance to note that due to the 
lack of self-interaction correction in common exchange- 
correlation functionals, the transferring (3d minority 
spin) electron will unphysically split between the two 
ions. Moreover, to calculate the energy gap, we need 
to accurately calculate the total energy when the minor- 
ity spin electron localizes at either reactant or product 
site at any given ionic configuration. 

In order to address these central problems, we first 
consider the simple case when oxidation states can be 
controlled trivially. This happens when two ions are in- 
finitely apart; the two ions can be studied in separate 
simulation cells and the oxidation states are controlled 
by simply changing the total number of electrons. For 
this special case, we performed runs using 

Fc (2+r)+ and 

F e ( 3 - r )+ (with r=0.0, 0.25, 0.50, 0.75 and 1.0), each sol- 
vated with 31 water molecules in the unit cell. We then 
carried out one Fe 2+ and one Fe 3+ run in the trajectory 
generated with Fe^ 2+r ^ + , and one Fe 2+ and one Fe 3+ run 
in the trajectory generated with Fe' 3_r ' )+ . (When calcu- 
lating energies of charged systems in periodic-boundary 



<© i «• «- 4§ & * 

7J « f ''I C 



•Q O c w 


| 1 


l 




2% 
7 

> v. ( - 




f • 4 


• 


S o „ 


■O " 





9 - c 


9 r 





O J & 0? $ Cb ^ 

tf &2M « 3?» 9 

r>" <5 c 7> p ' 



Product state Reactant state 

z ft J=E/p/?)duct/-E/reactant/ 

FIG. 2: Procedure used to calculate the diabatic free energy 
surfaces for electron transfer: The reaction coordinate at each 
time step is calculated from the energy difference between the 
product and reactant states in the ionic configuration pro- 
vided by the sampling run. The phase space is explored via 
the umbrella sampling parameter r, determining the oxidation 
state of the ions. 



conditions, the Coulomb interaction of charges with their 
periodic images should be removed |2*lj . In practice, 
these errors cancel out when calculating the energy gap, 
which is the difference in energy between systems with 
the same charge.) 

Fig. |31 shows the resulting diabatic surfaces; the final 
result is obtained by integrating [2^] 



F{e) 



J2F r (e) 9r (e) 

r 

E5r(e) 



(2) 



where F r (e) is the slope of the free energy curve in the 
different sections, each characterized by an umbrella- 
sampling parameter r, and the weighting factor gv-(e) is 
< 8{e — e(r)) > r . Each simulation lasted 5 ps after accu- 
rate thermalization. For this special case, the trajectories 
generated with Fe^ 2+r ' + and Fe^ 3_r ^ + are independent; 
since each of them provides n data points, there will be n 2 
energy gaps, providing high statistics and a very smooth 
free energy curve that fits accurately a parabola, with a 
coefficient of determination (i? 2 ) of 0.9996, and a reor- 
ganization energy of 1.77 eV. Note that at the tail ends 
of each sampling region the statistical accuracy becomes 
lower - this explains the slight deviations from a parabola 
seen in Fig. 

As mentioned earlier, the self- interaction errors of 
most exchange-correlation functionals result in a dra- 
matic qualitative failure in describing ions in different 
oxidation states when more than one ion is present. This 
failure can be exemplified by the case of two iron ions 
in the 2+ and 3+ oxidation state. When such a system 
is studied - e.g. using PBE-GGA - the HOMO electron 
will split between the two iron centers, to decrease its 
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FIG. 3: Diabatic free energy surface for ferrous-ferric elec- 
tron transfer in the special case when two ions are infinitely 
apart. The solid curve has been obtained from first-principles 
molecular dynamics. The dashed curves are mirror images, 
and correspond to a parabolic fit of the data. Different shades 
indicate portions of the diabatic surface sampled with r=0, 
0.25, 0.5, 0.75 and 1. 

own self-interaction. This behavior takes place irrespec- 
tive of the chemical environment of the two iron centers; 
we observe it for two isolated atoms, two hexa-aqua iron 
complexes, or two ions fully solvated. 

We will show in the following that this failure can be 
corrected by adding a penalty cost to ground states with 
non-integer occupation of the ion centers . This same 
approach is also used to calculate the Marcus energy gap, 
where for a given configuration we need to determine 
both the correct ground-state energy (with the transfer- 
ring electron in the reactant electronic configuration) and 
the first excited state (with the transferring electron in 
the product electronic configuration). We use and vali- 
date the following penalty functional 

pi ffi-fho x 2 

E{{^}\ - E[{^}}+J2 7= / exp(-— jjds, 

~ crjV27T J-oo to I 

(3) 

where fs^.\ is the largest eigenvalue of the minority-spin 
occupation matrix on ion I (calculated in this work by 
projecting the minority-spin Kohn-Sham orbitals on the 
3d orbitals of an isolated iron atom), and $ is its target 
value. To determine the optimal parameters, we sep- 
arately calculated the minority-spin occupation matrix 
for either a ferrous or ferric hexa-aqua ion embedded in 
a dielectric continuum (e=78) (24|. Then, we determine 
the parameters in the penalty functional so that the oc- 
cupation matrices of the ferrous or ferric clusters are ac- 
curately reproduced once the two are studied in the same 
unit cell (P 7 =0.54 eV, $=0.95 and ct/=0.01 on the fer- 



rous ion and P 7 =-0.54 eV, $=0.28 and o- 7 =0.01 on the 
ferric ion). We note that the target occupations for the 
minority-spin are not chosen exactly one or zero, since 
the orbital hybridization between the iron 3d orbitals and 
the lone pairs of the water molecules contributes to the 
projection onto atomic orbitals (an alternative would be 
to use projections onto the 3d orbitals or the maximally- 
localized Wannier functions of a solvated iron, instead 
of an isolated atom). When calculating the energy gap, 
the penalty functional contributions are taken away from 
the total energy; in any case, these effects are negligi- 
ble since these contributions cancel out when calculat- 
ing energy differences. Different constraints or penalties 
have been recently proposed for density-functional calcu- 
lations [23, |2g, 123, 123 ; we found our choice particularly 
robust, but several variations on the theme can be envi- 
sioned. 

A first, qualitative validation of this penalty functional 
is performed examining the charge density obtained by 
subtracting from a calculation with a ferrous and a fer- 
ric hexa-aqua ion in the same unit cell that of an iso- 
lated ferrous hexa-aqua ion, and that of a ferric hexa- 
aqua ion (a dielectric continuum surrounds the two clus- 
ters to remove long-range electrostatic interactions be- 
tween them). When the penalty functional is applied, 
the charge density around the ions reorganizes itself so 
that it produces a charge density that is the exact su- 
perposition of that obtained from the two independent 
calculations. 

We can make our validation quantitative by calcu- 
lating the energy gap for the system described, using 
two different penalty-functional calculations that impose 
to the HOMO electron to localize first on one, then 
on the other ion. This energy gap can also be calcu- 
lated exactly with PBE-GGA using the "4-point" ap- 
proach J2t|, provided that all long-range electrostatic 
interactions are screened out. The four calculations 
involve Fe 2+ in two Fe(H20)6 geometries (A and B), 
and Fe 3+ in the same geometries; the energy gap is 
[E A (Fe 2 +) + E B (Fe 3 +) - E B (Fe 2 +) - E A (Fe 3 +)]. We 
choose one configuration in which the hexa-aqua ions are 
fully relaxed, and three carved out from random steps 
in the molecular dynamics simulations. The energy gaps 
we obtained with the penalty functional are 0.632, 0.569, 
0.769 and 1.027 eV, in excellent agreement with the "4- 
point" values of 0.622, 0.542, 0.769 and 1.012 eV. It is 
worth mentioning that the energy gap is an excited-state 
property of the system, and thus in principle outside the 
scope of density-functional theory, which is a ground- 
state theory. However, since the charge densities of the 
HOMO and LUMO do not overlap, we can argue that all 
that is required is a description of the charge density that 
is locally correct (the excited state has an electron locally 
in equilibrium around an iron, oblivious of the other iron 
ion where it could sit more favorably). 

With these tools, we determined the diabatic free en- 
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FIG. 4: Diabatic free energy surface for ferrous-ferric electron 
transfer when the two ions are 5.5 A apart. Different color 
shades indicate portions of the diabatic surface sampled r=0, 
0.5 and 1. The right dashed curve is the parabolic fit of the 
data and the left dashed curve its mirror image. 



ergy surfaces for two iron ions separated by 5.5 A and 
solvated in 62 water molecules, in periodic boundary con- 
ditions. 5.5 A was suggested 0,|3fJ to be the optimal dis- 
tance for electron transfer. We show our results in Fig. 01 
together with a parabolic fit to the data. The reorganiza- 
tion energy (A) that we obtain is 2.0 eV jUl) m excellent 
agreement with the experimental value of 2.1 eV [llj . 
The energy barrier AG «0.49 eV, about a quarter of A, 
as expected. We also note that since the structural and 
electronic configurations of all microstates are available, 
reaction mechanism (e.g. inner vs. outer sphere transfer) 
can be analyzed in detail by restricting the analysis to all 
configurations that have a Marcus gap close to zero. Full 
first-principles accuracy also implies that bond-breaking 
and bond-forming reactions can be followed in detail, or 
that protons can be explicitly introduced to study the 
pH dependence of the reaction |32l ]. 

In conclusion, we have demonstrated how it is possible 
to obtain Marcus diabatic surfaces from first-principles 
molecular dynamics, where the entire system is treated 
quantum-mechanically, with the accuracy and predictive 
power that this approach entails. The case when two ions 
are at a finite distance requires special care in dealing 
with self-interaction errors and excited-state energies. In 
response to these challenges, we developed and validated 
a penalty functional that is able to control the oxidation 
states of ions, and that describes accurately both the elec- 
tronic ground state and the first excited state where the 
electron is transferred to the other ion. This approach 
can be successfully applied to a wide class of oxidation- 
reduction reactions, in solution (as it often happens in 
electrochemistry or biochemistry) or in the solid-state 
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